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DYNAMICAL ERROR ANALYSIS IN ORBIT DETERMINATION SYSTEMS 


ABSTRACT 

In this report, a method is described in which the overall error or 
uncertainty involved in spacecraft trajectory determinations is used to 
define a time dependent error bound for the coordinates of the spacecraft. 
The differences obtained by comparing tracking data with theory reflect 
errors or uncertainties in the Hamiltonian of the system, and in turn 
form a type of Canonical ensemble considered in statistical physics. Such 
considerations then allows one to derive a set of virtual forces that ac- 
count for the uncertainties in the Hamiltonian which give rise to the cal- 
?! culated differences from the true orbit. 
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DYNAMICAL ERROR ANALYSIS IN ORBIT DETERMINATION SYSTEMS 


I. INTRODUCTION 

There has been much interest recently concerning error analysis of space- 
craft trajectory systems. Error analysis can be defined as the ability to describe 
the effect of inherent uncertainties of an overall trajectory determination system 
on the computational accuracy of the system. There are several sources re- 
sponsible for the presence of such uncertainties, the most important being due 
to the inability to correctly describe the physics of the problem. Hence, the 
mathematical modelling of the forces involved, and the physical data that is used, 
provide at best only a good starting point for calculation of an orbit. In addition, 
the type and quality of the observational data introduce further complications. 

The usual method of handling such matters is to perform what is called an orbit 
improvement or differential correction over many revolutions of an orbit, by 
comparing computed with observed data of the spacecraft and ’correcting’ or 
updating the initial conditions of the differential equations of motion. Here then, 
the observed data obtained by the tracking systems is considered to be the true 
or correct orbit data. In reality, the accuracy of the observational data will 
depend upon the tracking system that is employed. Even after fitting the orbit, 
the post convergence residuals do not account for uncertainties in coordinates 
of the tracking sites. Most important, even though the constants of integration may 
be determined well, the accuracy of the calculated orbit will still depend upon 
the accuracy of the orbit generator (Reference 1). 

In this report, a method is described in which the overall error or uncertainty 
involved in trajectory calculations is used to define a time dependent error bound 
for the coordinates of the spacecraft. 


II. CONSIDERATIONS OF STATISTICAL PHYSICS 

Let a system be characterized by a set of generalized coordinates and their 
equations of motion. In the canonical form we have, 


d t 


_ 

- 1 * 


dq K _ 3H 

dt -oP K 


( 2 . 1 ) 
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where H is the Hamiltonian function. In the theory of gases, statistical thermo- 
dynamics introduces the notion of phase space in which the whole gas is repre- 
sented by a single point with coordinates p and q (Reference 2). Associated with 
any physical system, is an ensemble of points in phase space, each of which 
represents a different possible state for the system. This set of points is the 
Gibbs ensemble and each representative point moves along a curve according 
to the Hamiltonian equations. It is assumed that such a curve is uniquely de- 
termined by a complete set (q 01 , • • • » Poi * . • .) of initial conditions and there- 
fore these orbits can never intersect. In addition, the points of the Gibbs ensemble 
contained within a closed hypersurface can never pass through that hypersurface 
if the surface moves according to Equations (2.1). When special symmetries are 
present such as when energy is conserved, each representative point is forced 
to move on a hypersurface or ergodic surface. An important form for the density 
of representative points corresponding to thermodynamic states is the Canonical 
ensemble and is given by 


p (E) r. , (2.2) 

where 0 and Q are two constants characterizing the distribution. This ensemble 
admits an energy fluxuation in the system ( 0 - E ), corresponding to a set of non- 
interacting ergodic surfaces in phase space among which the phase point moves. 
From this, the expectation value of any observable A is then given by 


<A: 


A p d 


r, 


(2.3) 


where dr is a volume element in phase space. We shall return to this point in a t 

later section. 

In a recent article, (Reference 3) it was shown by similar statistical consid- 
erations, that by relating the angular momentum to a quantity defined as the 
virial parameter through a derived quantity called the Vector Density Function, 
it is shown that for any conserative system, the periodic variations of an orbit, 
and consequently the cross track error, will time average to zero. 

In this report, we shall attempt to describe the motion of the system point 
or phase point between ergodic surfaces by letting the energy fluxuation of the 
Canonical ensemble now describe an ’energy fluxuation' of our orbit system due 
to the many sources of error inherent in it, that is, since a correct Hamiltonian 
is not presently possible because of the many difficulties described above, the 
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system point will be said to occupy any position on any of the ’ergodie surfaces’. 
These ergodic surfaces then represent energy fluxuations arising from all un- 
certainties. We only know the position described by (2.3) above given as the most 
likely by the orbit generator. At this point, we try to make use of observational 
data which is compared with our corresponding calculated data to define uncer- 
tainty on a predicted basis for the orbit system, in terms of our ’Canonical 
ensemble. ’ 


III. FUNDAMENTAL EQUATIONS 

Corresponding to a path described by an artificial earth satellite, the system 
point will describe such a conic section in phase space governed by equations 
(2.1). ’Fluxuations’ in the motion of the system point and hence the computed 
orbit can then be described by the following time derivatives of the orbital ele- 
ments (Reference 4), 


d a 

TFT 


/ 3\ * * 

2 (1 - c 2 ) 1 2 1, R / o sin v t T' + T' o cos v] 


(1 - c 2 ) 1 2 [R' sin v + T' cos v 4 T' cos E] 

dt W 


dS 2 /s~(l - S) r W' COS I // 
d (ft a) 1 2 (1 - e 2 ) 1 2 

d ^3 r W' si n y; 

dt (ft aS) 1 ' 2 (1 - e 2 ) 1 2 

d . 

... [a R' VS (I • o 2 ) cos v 

dt c (ijaS) 1 ' 2 (1 - e 2 ) 1/2 

-2T' a /S(l - e co s E) s i n v 
- T' a c /S (1 - e cos E) cos v si n v 

r a e (1 - e cos E) W' yT - S sin y] 
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i (1 - e 2 ) 1 ' 2 
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(3.1) 


where E, v , and y?< are the eccentric anomaly, true anomaly, and argument of 
latitude respectively. 


The parameters R', T', and W' are defined as the components of total forces 
acting on the orbit in the radial, tangential and normal directions. In the spirit 
of the above discussion, we now consider that these force components are not 
real but virtual, that is, these forces are responsible for the 'energy fluxuations' 
or spread of the ergodic surfaces in phase space. In other words, the computed 
orbit has deviated from the true or observed orbit because of the presence of 
these virtual forces. It is our task to determine them, and once this is accom- 
lished, the change in time of the six orbital elements can be obtained from equa- 
tions (3.1), which in turn are then inserted into the Vinti orbit generator to 
produce corresponding changes in spacecraft coordinates, ±AX, ±AY and ±AZ. 
The plus or minus signs indicate a coordinate bound, or error bound in the 
ordinates as functions of time, since they derive from an error or uncertainty 
in the energy or Hamiltonian of the system. 

At this point, we now assume that these errors reflected as virtual forces 
are described in terms of changes of the eccentric anomaly with time, that is, 


/ 

R' R E 

X' “ x E (3*2) 

W' ~WE 

where R, T, and W are coefficients to be determined. Now Kepler's law states, 

M r n ( t + ) - E - e s i n E, ( 3 * 3 ) 

where 3 is the Vinti parameter related to the time of perigee passage. 

Differentiating, 
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1 2 a 3 2 


d E n 

fit (l-o cos K) ( 1 « c ons !'i 


(3.4) 


Using this, equations (3.2) become, 


2 F 


T E 


W E 


R ,: 1 

2 a - 3 2 

^ <Jr 

(T - 

(> cos F) 

fl (2 

T, 1 

2 fl -3 2 

d 2 q T 

(i - 

c cos E) 

d t 2 

W/d 

2 a"3 '2 

- 

(1 

o cos E) 

<i t 2 


(3.5) 


per unit mass, where the q R , q T , and q w are taken to be displacements along the 
radial, tangential, and normal directions to the orbit. Therefore, 


In order to determine q T , 


d 2 q T 

(1-0 COS F) 

d t 2 

a- 1 2 a“ 3 2 

r|2 c 'r 

(1 - o cos E) 

c\ t 2 

,1 '2 a ~3 2 

_*-4» rl 

<> 2 % 

(1 - e cos E) 

d t 2 

u 1 2 a" 3 2 


q R , and q w , we refer to Figures 3.1, 3.2, and 3.3. 


(3.6) 
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Figure 3.1 


Here, 


q z q R sin r- 

q x = q R cos - cos ^ 

q y ” q R COS i* sin ft 3 


where " is the geocentric latitude, and • is the Vinti parameter associated 
with the geocentric longitude of the node. 
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Here 



q z - q T s i n I 


q x - q T cos I cos 



( 3 . 8 ) 


q y = q T cos I sin 



where I is the inclination of the orbital plane. 
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Figure 3.3 


Here 


/ 


q z - q w cos i 

q x . n w sin I cos J3 
q y q w sin I sin $ 


(3.9) 
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Now summing all cf the individual q x , q y , and q ? and identifying them as x, y, and 
i inertial respectively, we have, 


'X; /COS COS ' , COS I COS('“ f 
\ / 3 \ 2 


/ I / 


/• 3 j, sin I cos /? 3 \ /q \ 

\ 


cos sin /y , cos I s i n + /? 3 ^, sin I sin 


<3.10) 


/ 

z/ \sin s in I, 


cos I 


or simply, 


x 


Mr \ 

y 

- X 


z 1 
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Inverting, 


(3.11) 


Mr) 


r | 

C, T 

I 

1! 

y 

1%/ 

1 

l z I 


(3.12) 


Since the inertial coordinates are known functions of the orbital elements (Reference 
5), we have that the coordinates q R> q T and q w are also. As a result, Equations 
(3.6) can then be expressed in the form, 


T ( 1 - p cos F) [h' 2 d-A da / d a 2 \ ^ d T 

a 1 2 a - 3 2 |^\ Hoj 2 y d t yj ^-2y 9 a 


d j!i aq T 

S,3|/ Ht + \.rit2/ 3/3 

•* 3 
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(l-o cos E) ff ' 2 a + ( d 2 a \ f q R 

”T 7 2“-3 t 2“' \7Ttv 





w = 


(1 - c cos E) 
/*-! 2 a -3 y 2 



• V 


d a 


(3.13) 


The right hand sides of equations (3.13) are functions of R, T, and W, so that we 
have, 


T = f (R, T, W) 

R = g (R, T, W) (3.14) 

W = h (R, T, W) 


a system of transcendental equations where the unknowns are located on both the 
left and right hand sides. By examining equations (3.13) it is seer* that each 
quantity is known in terms of R, T and W. For example, substituting equations 
(3.2) into (3.1) we have, 




e 2 )~i/2 


[Re sin v + T + T e cos v] 


dE 
d t 


d e 
d t 



p fi - e' 2 cos E 
1 - e cos E 



+ T 


cos E + T 


(cos E _ e) 
1 - e cos E 


dE 
d t 


dS 2 yS (1 - S) r W cos (v + p 2 ) dE 
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d t 


3 - n (1 - o cos R) 

YiT a S (1 - e 2 ) 


[ s i 


in v cos % * cos v si n 
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1., -a 

d t ( — 2 2 
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(it n) 1 2 r]t 
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,d t d t. 
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2 d t 

Differentiating again with respect to time we have, 


d 2 a 
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2 e I /x 1 ^ 2 a ~ 3 " 2 

(1 - c cos E) 3 la - e cos E) 

- T Vl - e 2 sin F] - sin E [ 

Vl - e 2 
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e R Vl - e 2 sin E 


x T + T e (cos E - 
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(3.15) 


11 



* 


/ 


d 2 e /j. 1 ' 2 1 - e 2 r 1 

d t 2 a (1 - e cos E)^ 1a 3/2 (1 - e cos E) 
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(1 - e cos E) 



(3.16) 


From equations (3.12) we have that, 



‘a -f. 




(3.17) 


where t = 1 6 denotes a, e, . . . , fi 3 . Here B X/3t , B y/B'C, and Bz/B-t 

are given in Reference 6, or if one wishes to bypass calculation of the right 
ascension, by differentiating equations (15), (16), and (17) of page 34 of Reference 
5, which give, 
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where i - 1, 2, 3, and b 3i = 0 for i = 1, 2, and one for i = 3. For i - 4, (5, 
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To complete the right hand side of equations (3.13), we now can obtain ,.2 x /,'. K 2 , 
'y/^ f L and «' 2 z/*i fc? by either differentiating the first partials given by ‘ * 
Reference 6, or by differentiating equations (3.17) through (3.19) again For 
example, from (3.19) we have, 
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and similarly for h 2 y/7)i 2 and c 2 z/di 2 . 
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and similarly for <d 2 y/'di 2 and d 2 * z/3 The method for computing the co- 
ordinate error bounds can now be defined, since the right hand side of equations 
(3.13) are known in terms of R, T, and W. 

IV. ALGORITHM 

The procedure for computing coordinate bounds is as follows: 

1. Enter into equations (3.15) an estimate for R,T, and W based on space- 
craft geometry and nominal orbit conditions, in units of Newtons -force. From 
past studies, a good initial estimate for a Vinti solution of a Relay II type orbit 
would be about 10 7 Newtons for R and W, and approximately 10 " 6 Newtons for 
T since the tangential forces and in-track position errors are usually larger 
(References 1). 

2. The corresponding A a, A e, . . . , A ,8 are used in the Vinti orbit generator 

to produce equivalent coordinate differences AX. , AYj^ , and AZ where i refers 

to the observation. 
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3. For any given tracking data such as radio direction cosine for example, 

a corresponding A L. and A M. can be determined, which also can be called com- 
puted residuals. 

4. If the ’true’ residuals obtained by comparing theory with each observation 
are denoted as JR, and R M , then the standard deviation of fit is given by 



n 

—— [A I 1 .) 2 + (A M^ 2 ] 
n - 6 / i 1 1 


i - l 


>2 

(4.1) 


where A A . and AM. are the differences between true and computed residuals, 
or (AL. - R l ) and (AM. - R M ) respectively. 

5. We wish cr to be arbitrarily close to zero. If our criterion is not met, 
we return to equations (3.13), insert our initial estimates for R, T, and W into 
the right hand sides and compute a new set of R, T, and W on the left. These 
new or iterated values are entered in the equations (3.15) again, and steps (2) 
through (4) are repeated. 

6. This process is continued until a self consistent solution is found for 
the system of equations (3.13) which allows o- to be smaller than some pre- 
assigned, arbitrary, positive number. It is these values of R, T, and W that 
are accepted and used in conjunction with the set of equations (3.15) together 
with the Vinti orbit generator to then produce a set of corresponding error 
bounds as functions of time, that are given as X(t) ±AX(t), Y(t) ±AY(t), and 
Z(t) ± A Z(t). 
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VI. CONCLUSIONS 

Equation (2.3) can now be considered as follows: Although an integrating 
anomaly equivalent to temperature is not immediately clear for this classical 
system, the most likely or expectation value for a coordinate will depend upon 
the form of the Hamiltonian in the ensemble. As a result, if the probability 
density p contains the Vinti Hamiltonian, then the most likely values of the 
spacecraft coordinates are those generated by the Vinti program. In addition, 
the size of the error volume (±AX) (±aY) (±AZ) or 2 ] A X| j A Y 1 1 AZ | that 
bounds these most likely values are also determined by the Hamiltonian. If all 
interactions were known exactly, then obviously, the size of the cube or error 
volume would shrink to a point, namely that determined by the Hamiltonian, and 
the most likely values would then be the exact values for all time. 

This same analysis, it is felt, may be applied to spacecraft trajectory sys- 
tems in the pre-launch phase. A set of virtual forces R, T, and W can be de- 
termined by selecting and applying tracking data in this program particular to 
certain classes of orbits. The error bounds or volume produced can then be 
classified as nominal. 

The interest at present is to interface this program with the existing Vinti 
orbit determination system and to apply it to several satellites of current 
interest. 
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